'''
Reference:
# 图像锐化
[1] https://opencv-python-tutroals.readthedocs.io/en/latest/py_tutorials/py_imgproc/py_gradients/py_gradients.html#gradients
[2] https://www.cnblogs.com/babycomeon/p/13206130.html
# 同态滤波
[3] https://www.cnblogs.com/phoenixdsg/p/8414277.html
【4】 https://blog.csdn.net/qq_34725005/article/details/82691487
'''
import numpy as np
import cv2
import matplotlib.pyplot as plt
import math

def smoothFilter(img):
	# 均值滤波
	blur = cv2.blur(img,(5,5))

	# 中值滤波
	median = cv2.medianBlur(img,5)

	plt.subplot(131),plt.imshow(img, cmap = 'gray')
	plt.title('Input Image'), plt.xticks([]), plt.yticks([])
	plt.subplot(132),plt.imshow(blur, cmap = 'gray')
	plt.title('blur Image'), plt.xticks([]), plt.yticks([])
	plt.subplot(133),plt.imshow(median, cmap = 'gray')
	plt.title('medianBlur Image'), plt.xticks([]), plt.yticks([])
	plt.show()

def guassSmoothFilter(img, d_0=50, d_1=80, n=1):
	# 巴特沃斯低通，指数低通和梯形低通
	rows, cols = img.shape
	crow,ccol = int(rows/2) , int(cols/2)

	center_x = np.floor(rows/2)
	center_y = np.floor(cols/2)
	D = np.zeros((rows, cols))
	# 巴特沃斯低通滤波器
	H = np.zeros(img.shape)
	# 指数低通滤波器
	H_exp = np.zeros(img.shape)
	# 梯形低通滤波器
	H_t = np.zeros(img.shape)

	for i in range(rows):
		for j in range(cols):
			D[i,j] = np.sqrt((i-center_x)**2 + (j-center_y)**2)
			H[i,j] = 1 / (1+np.power((D[i,j]/d_0), 2*n))

	for i in range(rows):
		for j in range(cols):
			# D[i,j] = np.sqrt((i-center_x)**2 + (j-center_y)**2)
			H_exp[i,j] = np.exp(-0.374*np.power(D[i,j]/d_0, n))

	for i in range(rows):
		for j in range(cols):
			if D[i,j] < d_0:
				H_t[i,j] = 1
			elif D[i,j] > d_1:
				H_t[i,j] = 0
			else:
				H_t[i,j] = (1/(d_0-d_1)) * (D[i,j]-d_1)

	f = np.fft.fft2(img)
	# 将低频成分搬移到图像中心
	dft_shift = np.fft.fftshift(f)
	print(dft_shift.shape)

	
	butter_fshift = H*dft_shift
	butter_ishift = np.fft.ifftshift(butter_fshift)
	butter_img_back = np.fft.ifft2(butter_ishift)
	butter_img_back = np.abs(butter_img_back)

	exp_fshift = H_exp*dft_shift
	exp_ishift = np.fft.ifftshift(exp_fshift)
	exp_img_back = np.fft.ifft2(exp_ishift)
	exp_img_back = np.abs(exp_img_back)

	t_fshift = H_t*dft_shift
	t_ishift = np.fft.ifftshift(t_fshift)
	t_img_back = np.fft.ifft2(t_ishift)
	t_img_back = np.abs(t_img_back)


	plt.subplot(221),plt.imshow(img, cmap = 'gray')
	plt.title('Input Image'), plt.xticks([]), plt.yticks([])
	plt.subplot(222),plt.imshow(butter_img_back, cmap = 'gray')
	plt.title('Butter Spectrum'), plt.xticks([]), plt.yticks([])
	plt.subplot(223),plt.imshow(exp_img_back, cmap = 'gray')
	plt.title('exp Spectrum'), plt.xticks([]), plt.yticks([])
	plt.subplot(224),plt.imshow(t_img_back, cmap = 'gray')
	plt.title('ti Spectrum'), plt.xticks([]), plt.yticks([])
	plt.show()


def sobelDeal(img):
	# Sobel 算子
	x = cv2.Sobel(img, cv2.CV_16S, 1, 0)
	y = cv2.Sobel(img, cv2.CV_16S, 0, 1)

	# 转 uint8 ,图像融合
	absX = cv2.convertScaleAbs(x)
	absY = cv2.convertScaleAbs(y)

	# addWeighted(src1, alpha, src2, beta, 0.0)
	Sobel_img = cv2.addWeighted(absX, 0.5, absY, 0.5, 0)


	# plt.subplot(121),plt.imshow(img, cmap = 'gray')
	# plt.title('Input Image'), plt.xticks([]), plt.yticks([])
	# plt.subplot(122),plt.imshow(Sobel_img, cmap = 'gray')
	# plt.title('Sobel Image'), plt.xticks([]), plt.yticks([])
	# plt.show()
	return Sobel_img


def prewittDeal(img):
	# Prewitt 算子
	kernelx = np.array([[1,1,1],[0,0,0],[-1,-1,-1]],dtype=int)
	kernely = np.array([[-1,0,1],[-1,0,1],[-1,0,1]],dtype=int)

	x = cv2.filter2D(img, cv2.CV_16S, kernelx)
	y = cv2.filter2D(img, cv2.CV_16S, kernely)

	# 转 uint8 ,图像融合
	absX = cv2.convertScaleAbs(x)
	absY = cv2.convertScaleAbs(y)
	Prewitt_img = cv2.addWeighted(absX, 0.5, absY, 0.5, 0)

	# plt.subplot(121),plt.imshow(img, cmap = 'gray')
	# plt.title('Input Image'), plt.xticks([]), plt.yticks([])
	# plt.subplot(122),plt.imshow(Prewitt_img, cmap = 'gray')
	# plt.title('Prewitt Image'), plt.xticks([]), plt.yticks([])
	# plt.show()
	return Prewitt_img


def LaplacianDeal(img):
	# Laplacian
	dst = cv2.Laplacian(img, cv2.CV_16S, ksize = 3)
	Laplacian_img = cv2.convertScaleAbs(dst)

	# plt.subplot(121),plt.imshow(img, cmap = 'gray')
	# plt.title('Input Image'), plt.xticks([]), plt.yticks([])
	# plt.subplot(122),plt.imshow(Laplacian_img, cmap = 'gray')
	# plt.title('Prewitt Image'), plt.xticks([]), plt.yticks([])
	# plt.show()
	return Laplacian_img


def HomomorphicFilter(img, gammaH=2, gammaL=0.25, c=1, d_0=60):
	img_f = img.astype('float32')
	rows, cols = img_f.shape
	log_img = np.log(img_f+1)
	fft_log_img = np.fft.fft2(log_img)
	center_x = np.floor(rows/2)
	center_y = np.floor(cols/2)
	D = np.zeros((rows, cols))
	H = np.zeros((rows, cols))
	for i in range(rows):
		for j in range(cols):
			D[i,j] = (i-center_x)**2 + (j-center_y)**2
			H[i,j] = (gammaH-gammaL)*(1-np.exp((-1)*c*(D[i,j]/(d_0**2))))+gammaL

	#H = H/(rows*cols)
	log_img_2 = np.fft.ifft2(H*fft_log_img)
	img_3 = np.real(np.exp(log_img_2))
	print(img_3.dtype)
	img_4 = np.zeros(img.shape)
	cv2.normalize(img_3, img_4, 0, 255, cv2.NORM_MINMAX)
	img_4 = img_4.astype('uint8')

	# img_3 = img_3.astype('uint8')
	# img_3 = img_3 / 255.0
	print(img_4[:20,:20])
	print(img_4.dtype)

	# dst = cv2.equalizeHist(img_4)
	
	# img_3.astype('uint8')

	plt.subplot(121),plt.imshow(img, cmap = 'gray')
	plt.title('Input Image'), plt.xticks([]), plt.yticks([])
	plt.subplot(122),plt.imshow(img_4, cmap = 'gray')
	plt.title('Prewitt Image'), plt.xticks([]), plt.yticks([])
	plt.show()



if __name__ == '__main__':
	img = cv2.imread("../images/Lenna.png", 0)
	homomorphic_img = cv2.imread("../images/coin.png", 0)
	salt_img = cv2.imread("../images/circuit-board-salt.tif", 0)
	guass_img = cv2.imread("../images/ckt-board-gauss.tif", 0)
	

	# 平滑滤波
	smoothFilter(salt_img)
	# 频域低通滤波
	guassSmoothFilter(guass_img)

	# 锐化
	Sobel_img = sobelDeal(img)
	Prewitt_img = prewittDeal(img)
	Laplacian_img = LaplacianDeal(img)

	plt.subplot(221),plt.imshow(img, cmap = 'gray')
	plt.title('Input Image'), plt.xticks([]), plt.yticks([])
	plt.subplot(222),plt.imshow(Sobel_img, cmap = 'gray')
	plt.title('Prewitt Image'), plt.xticks([]), plt.yticks([])
	plt.subplot(223),plt.imshow(Prewitt_img, cmap = 'gray')
	plt.title('Prewitt Image'), plt.xticks([]), plt.yticks([])
	plt.subplot(224),plt.imshow(Laplacian_img, cmap = 'gray')
	plt.title('Prewitt Image'), plt.xticks([]), plt.yticks([])
	plt.show()

	# 同态滤波
	HomomorphicFilter(homomorphic_img)
